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We theoretically investigate the ground state of trapped neutral fermions with population 
imbalance in the BCS-BEC crossover regime. On the basis of the single-channel Hamiltonian, 
we perform full numerical calculations of the Bogoliubov-de Gennes equation coupled with the 
regularized gap and number equations. The zero-temperature phase diagram in the crossover 
regime is presented, where the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) pairing state governs 
the weak-coupling BCS region of a resonance. It is found that the FFLO oscillation vanishes in 
the BEC side, in which the system under population imbalance turns into a phase separation 
(PS) between locally binding superfluid and fully polarized spin domains. We also demonstrate 
numerical calculations with a large particle number 0(10^), comparable to that observed in 
recent experiments. The resulting density proflle on a resonance yields the PS, which is in good 
agreement with the recent experiments, while the FFLO modulation exists in the pairing field. 
It is also proposed that the most favorable location for the detection of the FFLO oscillation is 
in the vicinity of the critical population imbalance in the weak coupling BCS regime, where the 
oscillation periodicity becomes much larger than the interparticle spacing. Finally, we analyze 
the radio-frequency (RF) spectroscopy in the imbalanced system. The clear difference in the 
RF spectroscopy between BCS and BEC sides reveals the structure of the pairing field and 
local "magnetization" . 

KEYWORDS: quantum atomic gas, FFLO state, BCS-BEC crossover, imbalanced superfluid, Bogoliubov- 
de Gennes equation, pliase diagram, radio-frequency spectroscopy 



1. Introduction 

There is increasing interest in the investigation of neu- 
, tral Fermi systems with mismatched Fermi surfaces.^ 

■ The robustness of superfluidity against the "paramag- 
, netic" depairing is a longstanding fundamental issue that 

has captured the attention of researchers in various fields, 
ranging from condensed matter to color superconductiv- 
ity in dense quark matter.^ Recently, superfluid (SF) 
phases under population imbalance have been realized 
in a trapped Fermi gas,'^"* accompanied with the ma- 

, nipulation of the s-wave scattering length a by Feshbach 

' resonance. 

For the achievement of superfluidity in neutral atom 

■ systems, the Feshbach resonance is a key factor. Applying 
, an external magnetic field enables one to control the rel- 
' ative energy between the two channels of the scattering 

process: the open (scattering) and closed (bound) chan- 
nels. Two fermions distributed in hyperfine spin states la- 
beled as a—'\ , I form the Cooper pair via the weak attrac- 
tive interaction in the system with a negative a. In con- 
trast, for the positive a, the energy of the closed channel, 
i.e., the binding energy, is characterized as Ei, = — 1/Ma'^ 
for mass M, which leads to the stable formation of tightly 
bound molecular bosons. The composite bosons with 
a long lifetime turn into the Bose-Einstein condensed 
phase below the critical temperature. Hence, the manip- 
ulation of interatomic interaction continuously changes 
the superfluidity from the fermionic Bardeen-Cooper- 
SchriefFer (BCS) type to molecular Bose-Einstein con- 
densation (BEC) through the unitary limit on resonance. 
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i.e., the BCS-BEC crossover.^ 

In actual experiments on neutral atoms, the total par- 
ticle number N is conserved. Applying a radio- frequency 
(RF) field, one can control the population difference be- 
tween two hyperflne spin states, called the population 
imbalance. 



where Ng. is the number of spin a species. The negligi- 
ble contribution of the dipole-dipole interaction leads to 
the conservation of population imbalance P throughout 
the typical experimental time scale. It is known that in 
the presence of population imbalance, i.e., P ^ 0, the 
uniform superfluid state cannot be thermodynamically 
stable. Various candidates for the pairing state when 
Pj^O situation have been proposed, including the Fulde- 
Ferrell-Larkin-Ovchinnikov (FFLO) modulated pairing 
state, ^^'^^ the BCS-normal phase separation (PS),^'^ the 
breached pairing or Sarma state, the deformed Fermi 
surface superfluid (DFS),^^' "'^^ and the p-wave pairing 
state. ^''■^^ These proposed pairing states are robust even 
in the presence of the imbalanced spin density, i.e., the 
"magnetized" or imbalanced superfluid. Studies on the 
thermodynamic stability of such exotic pairing states 
have a long history.^' 

Recently, the superfluid phase diagram for the ho- 
mogeneous system has been extended to the BCS-BEC 
crossover regime by a number of authors. Here, the 
phase diagram is constructed in a plane of the tempera- 
ture T and the population imbalance P. The PS appears 
in the BEC side. In the deeper BEC limit, the homoge- 
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neous BEC superfluid of the boson-fermion mixture is fa- 
vored. ^^''^^ The FFLO state becomes thermodynamically 
stable in a narrow window in the BCS side. These studies 
take into ac;count only one particular form of the FFLO 
pairing that has a single center-of-mass momentum vec- 
tor Q, i.e., A(r) = Aqc**^ ''. However, it should be empha- 
sized that for an arbitrary value of P, the oscillation of 
the stable FFLO phase can be described using the mul- 
tiple vectors Q and -Q, i.e., A(r) = Aq cos (Q • r).^^'^^ 
This effect has not been considered in previous works, 
except for Ref.^^ which considers higher harmonics of Q. 
The generalized FFLO phase may compete with the PS 
in the remaining area of the phase diagram. 

The presence of a trap potential, which is used to cap- 
ture atomic gas in actual experiments, may lead to a dif- 
ferent situation. The simplest way to take account of the 
trap is to employ a local density approximation (LDA), 
which is achieved by replacing the chemical potential 
with the local quantity including the trap potential. The 
LDA calculation""'""'* predicts that the PS state, i.e., the 
BCS state surrounded by the spin-polarized normal do- 
main, is favored under a realistic condition with a large 
number of particles and a fully three-dimensional trap. 
However, we should mention that the LDA eliminates 
the FFLO pairing state, which is one of the possible can- 
didates for the ground state, because of the lack of the 
gradual spatial variation of the pairing field. Taking into 
account the gradient effect, full numerical analysis has 
been performed by several authors in the weak-coupling 
BCS regime^^"^^ and at the unitary limit.^^"^^ They pre- 
dict the stability of the FFLO oscillation, which cannot 
be described in terms of a single Q. 

The aim of this paper is twofold. The first goal is to 
clarify the ground state of trapped fermions with popu- 
lation imbalance. Previously in our series of papers,^''' '''' 
we presented numerical results in the weak-coupling BCS 
regime. The current work covers a wider region, including 
the BEC side of a resonance in the plane of the popula- 
tion imbalance P and the dimensionless parameter kpa, 
where fcp is the Fermi wave number. To address such a 
problem, we start with the Bogoliubov-de Gennes (BdG) 
equation, which includes a contribution from the gra- 
dient effect of the pairing field and describes the physics 
in the atomic scale ~A;p^. 

The second goal of the present paper is to discuss how 
the quasi-particle structure in the BCS-BEC crossover 
regime under population imbalance affects the RF spec- 
troscopy, which has been experimentally performed by 
Schunck et al.^ 

Also, we investigate how a large particle number N = 
O(IO^) changes the FFLO oscillation. In this work, we 
fully solve the BdG equation coupled with the regularized 
gap equation and number equation, where the contribu- 
tions from the higher energy are supplemented by the 
LDA. This hybrid calculation enables us to demonstrate 
the stability of the FFLO modulation at the quantitative 
level, which is comparable to the results of recent ex- 
periments.^"* The numerical results indicate that FFLO 
modulation exists even in the vicinity of a Feshbach res- 
onance l/kptt'^ —0.5, which has been already observed 
experimentally. ^ 



This paper is organized as follows. In § 2, we derive the 
BdG equation coupled with the regularized gap equation 
on the basis of the single channel model. Also, we show 
the numerical results on the l/kpa dependence of ba- 
sic physical quantities in the BCS-BEC crossover regime 
without population imbalance. We present the ground- 
state structures of the imbalanced system in § 3, where 
the spatial profiles of the pairing field and densities in the 
strongly interacting system (l/A;F|a| < 1) are displayed. 
In addition, we shall present a quantum phase diagram 
in the l/kpa-P plane. In §4, we present the numerical 
results on the local density of states (LDOS) and the RF 
spectroscopy for the imbalanced superfluid in the BCS- 
BEC crossover regime. The final section is devoted to 
conclusions and discussion. In addition, we give supple- 
mentary information, e.g., the derivation of the thermal 
Green's function and the gap equation, in Appendices A 
and B. 



2. Theoretical 
Model 



Formulation: Single-Channel 



2.1 Bogoliubov-de Gennes equation 

Let us consider a Fermi gas distributed in two hy- 
perfine spin states {a =t,i). The Fermi system across 
a broad Feshbach resonance, which is realized in ^'Li or 
atoms, can be well described by the single-channel 
Hamiltonian: 



n 



Jdrjdr' 



^Vi(r)if(°V.(r)<5(r-r') 



-[/(r-r')^|(r)^|(r')^i(r')^T(r 



(2) 



with the creation and annihilation operators of fermions, 
ipl{r) and tjJa{r)- The single-particle Hamiltonian is given 

by 

Hi°\r) = -^\/' + V{r)-„^, (3) 

where atoms with mass M are trapped by a harmonic 
potential V(r). Throughout this paper, we set h=kB = 1- 
The interatomic interaction potential is U{r~r') and the 
chemical potential of two species is ^i^i=^± Sfi, where, 
without the loss of generality, wc set // f > i.e., the spin 
up (spin down) is the majority (minority) component. 

Following the procedure described in Appendix A, the 
Bogoliubov-de Gennes equation is given as 



/CT(r) A(r) 
A*(r) -/C|(r) 



, (4) 



where the diagonal element is obtained from /Co-(r, r') = 
5{r — r')JCa{r) in eq. (A-3b). This equation describes 
the quasi-particle state with eigenfunction [u^, v^] and 
eigenenergy Ei, under the pairing field A and Hartree 
potential gp^. Here, the interparticle interaction is char- 
acterized by the s-wave scattering, U {r — r') = gS{r — r') , 
whose "bare" coupling constant is g = Ana/M with 
s-wave scattering length a. Hereafter, the interaction 
strength of the system is characterized using the dimen- 
sionless form kpa with the Fermi wave number kp = 
\/2MEy- Ey is the Fermi energy in a noninteracting 
Fermi gas, whose definition is given in § 2.2. 
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The bare coupHng constant, however, provides two sin- 
gular contributions to the BdG equation: (i) the ultravio- 
let (UV) divergence of the pair potential A(r) and (ii) the 
divergence of the Hartree potential at the unitary Umit 
/cpci— >±oo. It is known^^"'''* that the UV divergence can 
be renormalized by replacing the bare coupling constant 
g with the effective constant ^(r) in the gap equation. 
The explicit form of the regularized gap equation^^' ^° is 
given as 



A(r) 



g{r)^u^{v)vl{Y)J^, 



(5) 



where the renormalized coupling constant ^(r) is given 

by 



1 

W) 



1 

= - + 

9 



1- 



^F(r) ^kc{r) 



(6) 



2fce(r) fcc(r) - fcF(r) 

Here, the Fermi distribution function is = f{E^) = 
l/(e^'^/"^ + 1). The summation in eq. (5) is carried out 
for all the eigenstates with both positive and negative 
eigenenergies, whose details are described in Appendix 
A. The above gap equation is now free from the energy 
cutoff Ec- For Ec^Ep, the expression for ^(r) coincides 
with that obtained from the two-body T-matrix in the 
absence of the medium."''^ Here, fcF(r) and fcc(r) arc the 
local wave vectors defined by the local Fermi and cutoff 
energies, respectively: 



Ef{t)=i,-V{v), 



MM 
2M 



(7a) 



(7b) 



Note that even if the gap equation with the bare cou- 
pling constant g is singular, the effective constant ^(r) 
yields a nonsingular negative value in an extensive re- 
gion, ranging from the unitary limit to the deep BEG 
limit. 

For the second divergent behavior (ii), it is known^^'^^ 
that the divergent term at the unitary limit can be renor- 
malized if we consider the many-body contributions be- 
yond the mean-field self-energy, where the system be- 
haves as a Fermi liquid with an effective mass. Hence, 
we remove the singularity by neglecting the Hartree 
term, that is, the diagonal elements in the BdG equation 
(eq. (4)) are replaced by the single-particle Hamiltonian 
in eq. (3): /C^(r) = iri°)(r). 

In summary, the renormalized coupling constant in 
eq. (6) and the neglect of the Hartree potential lead to 
the regularization of the BdG formalism. The BdG equa- 
tion (eq. (4)) is self-consistently coupled with the gap 
equation (eq. (5)) and the number equation 



(8) 



where the particle density in each spin state is obtained 
from Eqs. (A-4b) and (A-5) by 



PT(r)=X^K(r)|V., 



Pi(r)=El^-WI'(l- 



(9a) 



(9b) 



This formalism is now free from any divergence and 
provides a qualitative expression for strongly interact- 
ing Fermi systems in the BCS-BEC crossover regime. 
Note that this equation within the single-channel model 
gives equivalent results to those obtained from another 
mean-field theory based on the fermion-boson model in 
the case of a broad resonance.^'* Also, in the deep BEG 
limit (l/fepa— »-hoo), the BdG equation (eq. (4)) can be 
mapped to the Gross-Pitaevskii equation with a small pa- 
rameter A/\n\,^^ where A describes the wave function of 
the condensed molecular bosons. As we shall show later, 
the chemical potential at the unitary limit 1 / fcpd — » is 
estimated as ii/Ep = 1 + p with /3 = —0.4 on the basis 
of the current theory, which is comparable to the recent 
experimental result of /3 = — 0.54.^'^^ 

In performing the numerical calculation, we compute 
the gap equation (eq. (5)) by the following hybrid proce- 
dure: A(r) = ABdG(r) + ALDA(r)- The first term is com- 
posed of the contributions from low-energy eigenstates 
ained from th(^ exact diagonalization of 
the BdG equation. The quantity Alda with the higher- 
energy contribution above e'^'^'^^ < E^ < Ec is supple- 
mented by the LDA, whose explicit expression is given 
by 

dp A(r) 



ALDA(r) = g{r 



f 



(BdG) (27r)3 2E{p, r) 
x[f{E^ip,r)) + fiE^{p,r)) 



1], (10) 



with pP'^^^ = \J2Me'^'^'^\ Here, we set E^^^{p,y) = 
£(p,r) T Sii with £(p,r) = VKp, r)]^ + | A(r)|2, 
e(p, r) = /2m, — fi, and Ec = pl/2m — (i. Also, the 
high energy contribution to each spin density is expressed 
within the LDA as p = p''a^^^^ + pi^^^^ with 



If 



dp 



x/(ii;Ta(p,r))+ 1 



(BdG) (27r)^ 

e(p,r) 



1 + 



e(p 



Ml 
».r)j 



^(P. 

/(-£x,T(p,r)) 



(11) 



^(P,i-), 

Note that the spatial variation of the pairing field is 
mainly determined by the contributions from the eigen- 
states with energy close to the Fcirmi energy, while the 
eigenstates with the higher energy may be described 
within the semiclassical approximation, i.e., the LDA. 
This hybrid procedure has also been used in the numeri- 
cal analysis of the thermodynamic quantities at the uni- 
tary limit.^^ 

2.2 Calculated system 

We numerically solve the BdG equation (eq. (4)), 
which is self-consistcntly coupled with the gap equation 
(eq. (5)). The theory takes account of the trap potential 
and the mismatch of Fermi surfaces 5ijl= (/x-f — At|)/2. At 
each iteration step, the chemical potential /x is adjusted 
to fix the total particle number defined in eq. (8). In the 
current work, we consider a cylindrical symmetric sys- 
tem with trap potential V{y) = ^Mcj'^r'^ (r^ = .t^ -|- y^), 
and impose a periodic boundary condition with period- 
icity Z = 3d {d~^ = VMlj) along the ^-direction. Un- 
der such cylindrical symmetry, the quasi-particle wave 
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functions are written as M^(r) = Uy(r)e*'^'''^+*^^^ and 
Vi,{r) = w,^(r)e'*-''''^"'"^^^-' with quantum numbers along 
the azimuthal- and z-axis: qg = 0, ±1, ±2, • • • and qz — 
0,±27r/Z, ±47r/Z,- • •. 

The BdG matrix in eq. (4) is then transformed by 
spatial discretization into a banded matrix with respect 
to the radial axis, which can be solved using the Lanc- 
zos/Arnoldi algorithm implemented in the ARPACK li- 
braries.^^ Throughout this paper, we use total particle 
numbers of iV = 3, 000 and 150,000. The corresponding 
Fermi energies in A = are given as E-p — 32llj and 154a;, 
respectively, using the definition Ep/uj — (SOTrn^/lS)^/^ 
with nz = N/Z. Throughout this paper, we set Ec^'^'^^ — 
150w — 4:.7Ep for the case of 3,000 atoms and 

^(BdG) ^ 200w = 1.3£;f for A^ = 150,000 to limit com- 
putation time. However, the higher-energy contributions 
up to Ec — lOOOw are supplemented within the LDA, as 
shown in Eqs. (10) and (11). 

2.3 Calculations for balanced systems 

Before discussing the numerical results for superfluid 
states with population imbalance, let us present the basic 
properties of balanced superfluids, described within the 
above mentioned mean-field theory. First, in Fig. 1(a), 
we show the pairing amplitude Ao/Ep as a function 





Fig. 1. (a) Maximum values of pairing amplitude in balanced 
population (P = 0) at T = and N = 3, 000 as a function of 
1/kpa. The dashed line is the BCS form described in the text. 
The inset displays the chemical potential shift. The dash-dotted 
line corresponds to one- half of the binding energy Ei,/2Ep. (b) 
The order parameter of the molecular bosons 2|^bec('" = 0)P 
(circles) is compared with the total density p(r = 0) (dashed line). 
The definition for ^bec is given in the text. 



of l/kpa, where Aq is the maximum value of the pair- 
ing field at zero temperatures. It is found that in the 
weak coupling limit l/A:Fa<— 1, the pairing field can be 
asymptotically described by the standard BCS relation. 



Ao/Ep = 8e- 



-2 — 7T/kp\a\ 



, while the pairing in the opposite 



limit becomes a wave function of tightly bound molecu- 
lar bosons, i.e., the order parameter of BEC. In this BEC 
limit, it is known that the BdG equation in the single- 
channel model can be mapped into the Gross-Pitaevskii 
equation for molecular bosons®^ where the fermionic 
chemical potential becomes one-half of the binding en- 
ergy of the pairs Ei, / Ep = —2 / (kpa)'^ , shown in the inset 
of Fig. 1 (a) . The order parameter of the molecular bosons 

is expressed as 4'BEc(r) = ^^A{r), corresponding to 

the total fermionic density 2|\E'BEc(r)P —p{r). As shown 
in Fig. 1(b), this asymptotic behavior can be confirmed 
by the direct calculation using the BdG equation. 

The intermediate region nearby 1 / kpa = is smoothly 
connected from the BCS to BEC limits. Then, the pair- 
ing field A obtained from the single-channel model de- 
scribes the order parameter composed of the fermionic 
pairs and the wave function of the molecular BEC. The 
pairing amplitude at the unitary limit is Ao/Ep = 0.7, 
which is overestimated in comparison with Ao/£'f = 0.5 
in the strong-coupling theory.^^'^^ It is also found that 
the coherence length is saturated toward the length scale 
of the interatomic spacing ^ofcF = 2£'F/Ao = 0(l) as kpa 
approaches the BEC limit 1/kpa^oo. 

3. Ground States in the Imbalanced System at 
T = 

3.1 Superfluid states in weak-coupling BCS regime 

We now consider the imbalanced case for A^ = 3, 000 
fermions trapped by a cylindrical potential at T — 0. Fig- 
ure 2(a) shows the spatial profiles of the pairing field at 
P = Q (dashed line) and P = 0.48 (solid line) in the weak- 
coupling BCS side of a resonance for 1/kpa = —0.52. 
For this coupling constant, it is found that the quantum 
phase transition from the superfluid state to the nor- 
mal state is induced at the critical population imbalance 
Pc = 0.61. It can be seen from Fig. 2(a) that in the cen- 
tral region labeled (I), the population imbalance does not 
affect the pairing field. In contrast, the superfluid pair- 
ing field is quenched in the outside region labeled (III). 
The pairing field in the intermediate region (II) yields 
the spatial oscillation, i.e., the FFLO modulation, where 
the amplitude gradually decreases toward the edge of the 
cloud. 

These characteristics of A(r) are reflected by the 
density profiles of each spin component displayed in 
Fig. 2(b), where we define the local population differ- 
ence, called the local "magnetization" , as 



m{r) = p|(r) -px(r). 



(12) 



Spin states in region (I) attract each other, and the mag- 
netization is excluded in order to obtain the full conden- 
sation energy. In the intermediate region (II) the gap 
function changes its sign, allowing it to accommodate 
the excess majority species. This is indeed a characteris- 
tic of the FFLO state; The accumulation of excess ma- 
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Fig. 2. Spatial profiles of (a) the pairing field at P = (dashed 
line) and P = 0.48 (solid line) and (b) the corresponding density 
profiles, where the dashed, dotted, and solid lines denote the 
majority and minority densities and the local magnetization, re- 
spectively. The inset in (a) shows A(r) (solid line) and pair) 
(dashed and dotted lines) with a logarithmic scale. In (c), A(r), 
m(r), and P|(r) at P = 0.58 are displayed with solid, dotted, and 
dashed lines. All results are at T = and 1/A;pa = — 0.52. 



jority species at r = results from the topological struc- 
ture of A(r). The quasi-particles across the FFLO node 
undergo a 7r-phase shift of the pair potential, allowing 
them to form a mid gap state that is spatially bound 
there, called the Andreev bound state in more general 
contexts."^*' The energy of this state is situated in 
the middle of the energy gap, and the mismatch of the 
Fermi surface causes the difference in the occupation 
of this bound state, leading to the local magnetization 
around the nodes, as seen in Fig. 2(b). Here, a sufficient 





Fig. 3. (a) Spatial distributions of |A(r)| at l/fcpa = — 0.52 and 
T = for various values of P. (b) Corresponding local magneti- 
zation m(r) (solid lines) and the majority density P|(r) (dashed 
lines). The circles correspond to the radii of the pairing field, Rc 
defined as A(i?c)/Ao = 10"^. 



amount of the minority component remains for pairing 
with the majority component, i.e., the local magnetiza- 
tion is not fully polarized. These features are revealed by 
the bimodal distribution of the minority component (see 
around r/d ^ 3 in Fig. 2(b)). 

Figure 2(c) shows the gap and density profiles in the 
vicinity of the critical population imbalance P/Pc = 0.95. 
Here, the local magnetic moment r7i(r) governs the en- 
tire region of the system, where the "empty core" in the 
central region of the local magnetization vanishes. This 
leads to the quenching of the balanced BCS pairing even 
at r~0. As seen in Fig. 2(c), however, the superfluidity 
remains robust up to the edge of the minority component 
r~6.5(i by forming the FFLO pairing, accompanied with 
partially polarized spins. 

In Fig. 3 we display the P-dependence of the pairing 
field and local magnetization at r = and l/A:Fa = — 0.52. 
With increasing P, the oscillating region becomes wider 
towards the central region, while the radii of the pair- 
ing field Rc keep a constant value Rc/d ~ 7 up to the 
critical population imbalance of Pc = 0.61. This FFLO 
pairing state is not describable with the LDA where the 
superfluidity is localized in the central region (I) and 
the regions (II) and (HI) are regarded as being in the 
normal state, i.e., the BCS-normal PS state. It is 
found that this oscillating pairing state becomes robust 
up to the critical population imbalance Pc = 0.6 at r = 
and l/fcpa = —0.52; The equal population P = is the 
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only stable situation for the nonoscillating BCS state. Be- 
yond Pc, the superfluid state becomes the normal state 
through a second-order phase transition. 

In the outside region of Rc, where the gap almost 
vanishes, the complete spin-polarized state is attained. 
It should be emphasized again that the central region 
r ^ 0, in which the magnetization is completely excluded, 
catches a clear signature of the "balanced" superfluidity, 
while the surrounding area with the partially polarized 
spins also keeps its superfluidity composed of the "im- 
balanced" FFLO pairing. 

3.2 Strong- coupling region: from unitary limit to BEC 
regime 

Let us now consider the strong-coupling BEC region. 
The pairing field in the BEC side of the resonance is 
completely different from that in the BCS side. Figure 
4 shows the spatial profiles of A(r), Pa{r), and m{r) at 
l/fcFa = 0.52, P = 0.73, and r = 0. The superfluid pairing 
state is still robust in the central region of the system, 
while the outside region turns to the normal state. The 
intermediate region (II) smoothly connects the superfluid 
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Fig. 4. Spatial profiles of (a) the pairing field at P = (dashed 
line) and P = 0.73 (solid line) and (b) the corresponding density 
profiles, where the dashed, dotted, and solid lines denote the 
majority and minority densities and the local magnetization, re- 
spectively. The inset in (a) shows A(r) (solid line) and pair) 
(dashed and dotted lines) with a logarithmic scale. All results 
are for the BEC side l/%a = -|-0.52 at r = 0. 



core (I) and the normal state (III) without any spatial 
oscillation of the pairing. In contrast with the BCS side 
shown in Fig. 2, it is seen from the inset of Fig. 4(a) that 
the intensity of the pairing field exponentially decays to- 
ward the edge of the minority component. As seen in 
Fig. 4(b), the magnetization is completely excluded from 
the central region (I), while the spins in the outer region 
(III) are fully polarized. In the strong-coupling regime, 
the intermediate region (II) plays the role of the domain 
wall between the superfluid core (I) and the fully po- 
larized normal state (III). This phase-separated profile 
is commonly seen in other theoretical calculations based 
on the LDA,'^""'^^'^^''^^ except for the proximity effect of 
the pairing field in the vicinity of the domain wall. 

Figure 5 displays the spatial distributions of the pair- 
ing field and the local magnetization as a function of 
P, where the circles denote the edge of the pairing field 
i?c defined as A(i?c)/Ao = IQ-^. in the BEC side, the 
critical population imbalance is uniquely determined as 
Pc = 1 at zero temperatures. The robustness of the su- 
perfluidity is supported by the fact that it is realized by 
the formation of local pairs, i.e., molecules. 

Also, for the P-dependence of the ground state shown 
in Fig. 5, two clear differences exist between the BEC 
and BCS regimes; First, in the case of l/fcFa = 0.52, the 
radii of the condensation area Pc gradually shrinks as P 
approaches Pc = 1, while Pc in the BCS regime has an 
almost fixed value up to Pc. Second, a domain composed 
of "fully" polarized spins grows around the empty core 
for all values of P in the BEC side. In contrast, as seen in 




Fig. 5. (a) Spatial distributions of |A(r)l at the BEC side 
l/fcFa = 0.52 and T = for various values of P. (b) Corresponding 
local magnetization m(r) (solid lines) and the majority density 
Pl(r) (dashed lines). The circles correspond to the edge of the 
pairing field Rc defined as A(_Rc)/Ao = 10"'^. 
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Fig. 3, the outside area in the BCS side has a sufficiently 
large region (II), which is composed of "partially" po- 
larized spins and the minority spins to take part in the 
superfluid pairing with the majority spins. 

The pairing field and density profiles at the unitary 
limit are displayed in Fig. 6 where the FFLO oscilla- 
tion remains in the outside area of the "core" , which has 
the balanced spin density. This modulation survives as 
a proximity effect between the balanced superfluid and 
the fully polarized domains, analogous to superconduc- 
tor /ferromagnet interfaces.''^ Note that the periodicity 
of the FFLO oscillation may be scaled with the coher- 
ence length ^0- As shown in the inset of Fig. 6, the FFLO 
oscillation is completely periodic with periodicity L^d, 
which is comparable to ^o = 2.5A;p^=0.31d, i.e., L~3^o- 
In the BCS side for kpa = -0.52, the FFLO modula- 
tion region (II), which becomes wider toward the edge 
of the cloud, has longer periodicity, L~2rf = 2.7^0) e-g-> 
see Fig. 2, where = S/cp^ = 0.75rf. In contrast, the 
coherence length saturates at the interparticle spacing. 



Fig. 6(b) is almost unchanged compared with that in 
the BEG side (See Fig. 4), indicating a phase separation. 
This is different from the density profile in the BCS side 
as shown in Fig. 2(b). 

3.3 Large-N system 

We now turn to the situation for the realistic particle 
number TV = 150, 000. The spatial profiles of the pairing 
field and densities at l/fcpa^— 1-2 and —0.7, and in the 
vicinity of the unitary limit 1/A;fci = — 0.14 are displayed 
in Figs. 7 and 8, respectively. The tendencies due to the 
FFLO pairing, which have been already seen in the sys- 
tem with TV = 3, 000 atoms, are commonly reproduced 
even in the large- A'^ system. For instance, the outside re- 
gion r > Sd of Fig. 7 shows that the pairing field yields 
the FFLO oscillation, whose indirect signature at the 
macroscopic level is the partially polarized spin density. 
Note that recent experiments^ have been performed over 
the wide kpa range, l/Zcplal < 0.5, in which the FFLO 
modulation survives. The periodicity L of the oscillation 
in Fig. 8(a) has a larger length scale than the interpar- 



~ ' BEG side, in which the oscillation sud- 

denly vanishes. The resulting density profile shown in tide spacing, L ~ 4d = 70A:p\ which is comparable to 

the coherence length ^q: L = 3.3^o with S,o = 21k^^ and 



--OMd. 




The spatial profiles in the strong-coupling regime are 
shown in Fig. 8. The local magnetization is in good 



Fig. 6. Spatial profiles of (a) the pairing field ami (b) the cor- 
responding density profiles at P = 0.77 at the unitary limit 
l/fepa = 0. In (b), the dashed, dotted, and solid lines denote 
the majority and minority densities and the local magnetization, 
respectively. The inset in (a) shows A(r) (solid line) and Pa(r) 
(dashed and dotted lines) with a logarithmic scale. All results 
are at T = 0. 




Fig. 7. Spatial profiles of the pairing field (solid line), the ma- 
jority spin density (dotted line), and the local magnetization 
(dashed line) in the system with Ar = 150, 000 at (a) l/feFO = — 1.2 
and P = 0.13 and (b) l/fcFa=-0.7 and P = 0.12. All results are 
at T = 0. 
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Fig. 8. Spatial profiles of (a) the pairing field and (b) the density 
in the system with Ar = 150, 000 at l/A:Fa = -0.14 and T = 0. The 
dotted, dashed, and solid lines in (a) show the pairing fields at 
P = 0, 0.06, and 0.29, respectively, and those in (b) denote the 
minority and majority densities, and the local magnetization, re- 
spectively. The dot-dashed line in (b) corresponds to the density 
of the up spins at P = 0. In the inset of (a), |A(r)| at P = (dot- 
ted line) and 0.29 (solid line) are displayed with a logarithmic 
scale. 



agreement with the results of an experiment in which 
the magnetism was observed by phase-contrast imaging 
and 3D image reconstruction.^ While the local magne- 
tization yields a PS-like profile, FFLO modulation oc- 
curs in the pairing field. In the strong-coupling regime 
l/kpa^ —0.14, the periodicity of the FFLO oscillation 
becomes shorter, which is comparable to the interparti- 
cle spacing L ~ 0.6(i~ 10.5A;p^, and its intensity expo- 
nentially decays as seen in the inset of Fig. 8(a). With 
^o = 3.7fcp^, the oscillation period is scaled as L^2.8^o- 
In summary, throughout the extensive range of the in- 
teraction fcpfl < 0, the oscillation periodicity L is well 
scaled with the coherence length as L = a^o- The co- 
efficient a is around 3. We also find that this result is 
insensitive to the total particle number, i.e., the Fermi 
wavelength k^^. Surprisingly, it is found that the period 
is almost unchanged with an increase in the population 
imbalance. One example of this is displayed in Fig. 3(a), 
where internode spacing is almost fixed for the entire 
range of P. Note that since the FFLO oscillation peri- 
odicity in the absence of the trap potential is sensitive 
to the population imbalance or alternatively to the mis- 



match of the Fermi surface,"^*' the constancy of L may 
be peculiar to the finite trap system. 

3.4 Quantum phase diagram 

Let us now summarize the ground state of the imbal- 
anced Fermi system in the BCS-BEC crossover regime 
by constructing the phase diagram at zero temperatures. 
Figure 9 shows the phase diagram in the plane of the di- 
mensionless coupling constant kpa versus the population 
imbalance P for the system with TV = 3, 000 atoms. It is 
important to mention that the critical population imbal- 
ance in the weak-coupling limit 1 /fcpa < — 1 exponentially 
depends on l/fcpa, i.e., P^oce^'^^'^''^^'^^ which results in 
the linear relationship with the intensity of the pairing 
field, i.e., Pc oc In a previous work,^" we found that 
Pc = 1.9^. In Fig. 9, the corresponding exponential line 
is depicted with 1.9Ao/i?F where Aq is the maximum gap 
A(r = 0) at r = and P = 0. It is seen that this tendency 
is also confirmed in the current phase diagram obtained 
from the crossover theory. The superfluid phase below 
Pc at fcFa<0 yields the spatial oscillation of the pairing, 
i.e., the FFLO state. The formation of the FFLO pairing 
pushes up the phase boundary relative to the Pauli limit 
(see also Fig. 13). 

The behavior of the Pc curve in the BEC side is in con- 
trast with that in the BCS side, where the phase bound- 
ary is uniquely determined as Pc = 1. This results from 
the fact that the superfluidity survives by locally forming 
a molecular-like pairing with the corresponding amount 
of majority spins. The SF phase at 1 /kpa > corresponds 
to the PS state without any oscillation of the pairing. 
The distinct phase boundary between the FFLO and PS 
states cannot be defined because the FFLO state contin- 
uously turns into the PS state via the proximity effect 
in the BCS/polarized-normal domain interfaces, as has 
been discussed above. This is peculiar to the finite trap 
system. 

It has been proposed'^^' that as l/fcpa reaches the 
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Fig. 9. (Color online) Quantum phase diagram in l/fcpa-P plane. 
The points denote the estimated phase boundary between the 
SF and normal (N) phases. The pairing field in the superfluid 
phase exhibits the spatial oscillation (nonoscillation) for negative 
(positive) l/fcpa. The dashed-dotted line is the extrapolation 
from the results of the weak-coupling limit, Pc = 1.9Ao/-Ef,^'' 
where Aq is the maximum intensity of A(r) at T = and P = 0. 
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deep BEC limit l/fcpa ^ 1, the PS state becomes a 
"homogeneous" imbalanced superfluid without a phase- 
separated domain, that is, a mixed system of bosons and 
spinless fermions is attained. 

4. Quasi-particle Structure and Radio- 
Prequency Spectroscopy 

4-1 Local density of states 

The LDOS for each spin component is given by the 
definition 

K{r,E) = --Qg^{rr,iu}n^E + ir]), (13) 

TT 

where 5|(rr',za;„) = CJii(rr', iw„) and 5|(rr',zu;„) = 
— ^^22(r'r, — iw„). Using the thermal Green's function Qij 
described in Appendix B, one can read the LDOS for the 
spin-up state, 



E 



u,{r)\^6{E-E,), 



and for the spin-down state. 



^f^{r,E) = J2\Mr)\''S{E + E,). 



(14b) 



The LDOS for majority spin states in the BCS side 
are displayed in the top row of Fig. 10. In the balanced 

case, the low-lying excitations are bound in the surface 
of the cloud (r > 4d). It is known^^"^^ that the quasi- 
particles around the surface experience the effective po- 
tential consisting of the trap potential V{r) cx and 
the pair potential A(r), where the latter is a monotoni- 
cally decreasing function in the surface region. This sit- 
uation may be reduced to a problem on quasi-particles 
confined in the quantum well like potential A(r) -I- V{r). 
Hence, the eigenenergies close to the Fermi level are dis- 
cretizcd using the trap unit, i.e., a finite small energy 
gap of size O{oj) exists. For finite P, however, the FFLO 
oscillation induces a mid gap mode with zero energy, and 
the surface excitation gap vanishes as P increases. The 
quasi-particles around the FFLO node behave as gapless 
normal particles, but the superfluidity survives. 

Note that the LDOS profile for the minority species is 
almost the same as that for the majority spins, except for 
the shift of the Fermi level, which is shifted downward 

The bottom row in Fig. 10 shows the LDOS in the 
BEC side. At l/fcpa = 0.52, the maximum value of the 
energy gap is comparable to the Fermi energy |Ao| ^E'f- 
In contrast, the quasi-particle structure in the surface re- 
gion is different from that in the BCS side. In the deeper 
BEC limit, the low- lying excitation can be characterized 
by the binding energy, that is, E^ = \J /i^ + ~ |^|. 
Also, the quasi-particle states in energy bands lower than 
— A(r) — 6n are no longer the eigenstates of the harmonic 
oscillator, in contrast with those in the weak-coupling re- 
gion, l/fcFa = — 0.52. With increasing P, as shown in the 
bottom row of Fig. 10, the particles, for instance, for 
rjd > 4 at P = 0.56, locally dissociate from the pairing 
state, which becomes the normal state. 



4-2 Basic formalism for radio-frequency spectroscopy 

The spatially averaged quantity of the LDOS is observ- 
able in the RF spectroscopy, which has been recently de- 
veloped by several experimental groups, in an equal mix- 
ture^^"^^ and an imbalanced system.* Theoretical stud- 
ies have been carried out by a number of authors based 
on linear response and nonlinear response the- 
ories. ^^'^^ Here, following their previous works based on 
linear response theory, we describe the formulation for 
RF spectroscopy in an imbalanced system. 

We consider another internal state |e) in addition to 
the two hypcrfinc states \a =t)i)) which form a pair- 
ing via an effectively attractive interaction. The state |e) 
can be described by the field operators V''e(r) f^^nd %i'l{r), 
which obey the standard fermionic commutation relation 
and are commutative with those of other internal states. 



On the basis of several works 



73, 74, 80-84, 86 



extend the 



(14a) original Hamiltonian describing the state \a) to 



Ti. = Ti-MF + We + Ti-T + 



with 



Afa= dripl{r)ijja{r), a = <7,e. 



(15) 



(16) 



The Hamiltonian includes the following contributions. 
First, the Hamiltonian for atoms in the e state is given 

by 

He = / dr^l(r) Inf^ + J2 geaPair)] V'e(r), (17) 



with Tie = — + V{r) — /i,. . Second, the Hamiltonian de- 
scribing the "tunneling current" between internal states 
is introduced as 

n!f^ = Jdr [l](r)V|(r)Va(r) + h.c] . (18) 

The detuning frequency Wdet expresses the difference be- 
tween the internal energy level difference and the fre- 
quency of the applied laser. 

We calculate the current from the pairing state \a), 
corresponding to the rate of change of the population of 
the e-state: I{t) = {J\fe{t)). Here we divide the Hamil- 
tonian in eq. (15) into two parts: (i) the diagonal part 
Ho =Hmf +'He + ^ [Ea - ^e] and (ii) the pertur- 
bation Hamiltonian Ht- From linear response theory, 
the tunneling current is given as 



/(<^)(t) 



(19) 



where the expression of the time-dependent quantity is 
given by ©(t) = e*^"*Oe~*-^o* with the canonical Hamil- 
tonian Hq = Hq + HeNe + Ecr \J^oMa- In particular, the 
tunneling Hamiltonian defined in eq. (18) is transformed 
to 

W^^^W = / c;rfi(r) [e-*'^*Vl (r, i)V'a(r, t) + h.c] , (20) 
using the field operators in the Heisenberg represen- 



tation, Ve(r,f) 



V^e(r)e- 



and V<7(r,i) 



giMMFt^^j-j-^g iWMFt Hereafter, we use the notation cDe 
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Fig. 10. (Color online) LDOS for majority species for various values of P. The top row is in the BCS side (l/k-pa, = —0.52) and the 
bottom row is in the BEC side (l/fcFa = +0.52). The origin of the vertical axis corresponds to the energy equal to the chemical potential 
/i, and the shift of the gap center from the origin is characterized by the mismatch The solid lines denote the local energy gap 
defined by ±|A(r-)|. 



Wdet + Mo- ~ A'e for convenience. Also, the number op- 
erator for the e-state obeys the Heisenbcrg equation, 
iAfe{t) = [N'e{t),H], which leads to the following expres- 
sion for the Heisenberg representation, 

K{t) = -i j drn{r) [e'^^'Vl (r, OV'o-Cr, i) - h.c] . (21) 

By substituting Eqs. (20) and (21) into the Kubo for- 
mula in eq. (19), one can see that the total current is 
composed of two contributions, a single-particle tunnel- 
ing current /^'^■' and a Josephson tunneling current /j'^^ , 

as I^'^-' = /^'^■' + I'f^- Since we are interested in single- 
particle tunneling, the Josephson current is neglected 
here. The resulting single-particle tunneling current is 
obtained from the analytic continuation of the quantity 
expressed as the product of the Green's functions for the 
a- and e-states:*^ 

I^^\t) 23 y" y" W^(rr',iw„ ^ w -f iTj)drdr' , (22a) 

where the Matsubara frequency is introduced as ujn = 
7r(2n + 1)P and 



The Green's function for the e-state is given as 



^^4(rr^^u;„)=r'^^*(r)^^(r')^ 

X (rr' , iw„ - iw' J (rr' , ZLJ^ ) . 



(22b) 



(22c) 



C/e(rr',iLj„) = 



0c(r)0J(r) 



C 



'J71 — ec 



(23) 



where the eigenfunction and energy, 0^ and e^, are ob- 
tained from the Schrodinger equation for atoms in the 

e-state, [He"'' +J2 gecrP<T](f>ci^) = ^C'^c(^)- The expressions 
of the Matsubara Green's function for the pairing state 
Qii and Q22 are shown in Appendix B. 

To this end, one can find single-particle tunneling cur- 
rents for the following distinguishable processes: for the 
tunneling from the majority species \a =t) to the e-state, 



= 27r^ f O(r)M^(r)0*(r)di 



x[U^fiec)]Siu + E,-ec), (24a) 

and for the tunneling from the minority species |(t =J,) 
to the e-state, 



Afu- f{ecm^-E,-ec). (24b) 

Note that similarly to the BdG formalism in § 2, the 
summation in Eqs. (24a) and (24b) is carried out for 
all eigenstates with both positive and negative energies 
because of the breaking of the time-reversal symmetry. 
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In performing the numerical calculation, the i5-function 
in eq. (24) is replaced with the Lorentzian function 
d{z) Tniz) = (??/2)V[02 _^ (r;/2)2], where the resolu- 
tion of the spectrum 77 is set as 77= l.Ow throughout this 
paper. We consider the situation when the interaction 
between the pairing state and the e-state is negligible, 
fleer = 0. We also focus on the transition from the pair- 
ing state |(t) to the excited state |e), i.e., the positive 
detuning Wdet > and Ig'^ > 0. 

4-3 Numerical results 

It is important to mention that the tunneling current 
in the homogeneous pairing field A(r) = A at P = is 
expressed as 



pqi 



Me + Wdet) 



-/(e - t-ie)Wa{ci, e - HaWe{p,C + t^dot - Me), (25) 

with fipq = / dre*P ""J7(r)e-*i ''. The density of states for 
the pairing state is 

K (q, z) = 27r [ul6{z - E^) + v^Siz + E^)] , (26) 

and that for the e-state is Me{z) = 2-k5{z — Cp). Here, 
(uq,Uq) and ii'q are the solutions of the BdG equation 
in the homogeneous system at P = 0, and €p=p'^/2m — 
Me is the eigenenergy of a free particle in the e-state. 
We are interested in the situation of positive detuning 
and current, LO^et > and Ig'^ > 0, which expresses the 
fraction loss of the pairing tr-state. First, one obtains the 
following simple result in the case of A = 0, 

4"^ ('^det) = 27rO' {N^ - Ne) S{0Jdet), (27) 

which leads to a single-peak structure at Wdct = when 
the e-state is initially not occupied (/Xe = 0). The inten- 
sity of the peak gradually decreases with increasing n^, 
and in the case of the equal chemical potential 11^ = l^-a, 
the tunneling current is not responsible to any detuning 
frequencies, i.e., Ig'^ =0. 

For the pairing state A 7^ 0, after evaluating the inte- 
gral in eq. (25) over energy and momenta, one obtains^""^ 



det 



Wdet 
2 



+ 2/i„ 



A2 



'"det 



(28) 



where we introduce the density of states in the ideal 

Fermi gas, AfaiE) = -^Ve. Also, we set 5fl = fia — Me- 
For positive detuning Wdet > and a system with posi- 
tive chemical potential Mct >0, it is found that the func- 
tion Ig'\u)det) becomes monotonically decreasing in the 
range uJdet/Ep > A/Ef when the e-state is initially oc- 
cupied, fie = fJ-a- In contrast, in the case of Me = 0, the 
resonant detuning at which J^''' has the maximum value 
is shifted to Udct/Ep = |(.^)2. In the deep BEG limit 

where Mo- < and A/l/x^l "Cl, Ig'^ exhibits different be- 
havior from that in the BCS side, which is insensitive 
to /ig. Then, the resonant detuning is situated around 
Wdet/-E'F~2|/XCT|/-£'F~ l-Efel/i^F- This energy corresponds 
to the dissociation of molecular bosons and is uniquely 
characterized by the dimensionless parameter l/(fcFCi)^- 



In Fig. 11, we display the numerical results of the frac- 
tion loss of the minority species I^P in imbalanced sys- 
tems in the presence of a harmonic trap, where (a) and 
(b) correspond to the BCS and BEC sides, respectively. 

Hereafter, we consider the situation that the e-state is 
initially empty, Me = 0. Also, we set r2(r) = ri. In the BCS 
side, as shown in Fig. 11(a), the resonant detuning for low 
values of P is situated around o^dctZ-E-p — 0.1, which is re- 
lated to the dissociation energy of the fermionic pairing 
described above, uidet/Ep ^ (Aq/Ef)^, with the maxi- 
mum gap of Aq = 0.34i?F = llw. As P increases, the 
peak position approaches the zero detuning. Note that 
for high values of P, additional fraction loss occurs at 
the zero detuning and the resulting spectrum profile has 
a double-peak structure. This reveals the fact that the 
system under high imbalance is in the pairing state with 
partially polarized spins, which is indirect evidence of the 
FFLO state. As seen in Fig. 10, the mid gap state ap- 
pears in the spacing between the small energy gap near 
the surface {r/d ~ 6) when P 7^ 0. The presence of the 
mid gap state increases the intensity of the fraction loss 
at Wdct = 0. 

It is seen from Fig. 11(b) that the fraction loss in the 
BEC side yields a spectrum distinctive from that in the 
BCS side. There are two differences, (i) The detuning 
at which the fermionic/molecular pairs are dissociated 
is unchanged for increasing values of P. Only the inten- 




Fig. 11. RF spectroscopy of the minority component Ig'-' at (a) 
1/A;fo = -0.52 and (b) l/fcFa = 0.52 -with AT = 3, 000 atoms. All 
the results are at T = 0. 
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Fig. 12. Fraction loss of (a) majority and (b) minority species at 
l/fepa = —0.14 with AT = 150, 000 atoms. All the results are at 

r=o. 



sity becomes weak, (ii) With increasing P, no additional 
peak appears around the zero detuning. This is because 
all the minority spins in the imbalanced situation form 
the "pairs" with the corresponding amount of the major- 
ity species in the local region, and the spins in the normal 
state are fully polarized. The fraction loss for the major- 
ity species is unchanged in the extensive region from the 
BCS limit to the BEC limit, where the spectrum always 
yields the double-peak structure with one peak situated 
at LOdet = having large intensity and the other depend- 
ing on the dissociation energy of fermionic or molecular 
pairs at Wdet — ^o- 

Finally, the RF spectroscopy for the system with a 
realistic number of particles = 150,000 in the vicin- 
ity of the unitary limit l/k^a = —0.14 is presented in 
Fig. 12. The qualitative behavior is unchanged from the 
case of the small particle number displayed in Fig. 11, 
that is, the fraction loss of the minority species occurs 
at zero detuning in addition to ti;aet/-E'F~ |(^)^ = 0.17, 
corresponding to the dissociation energy of the pairing 
state at r = 0, Ao = 0.53i?F = 82a;. The corresponding 
pairing field is shown in Fig. 8, where the FFLO mod- 
ulation appears in the vicinity of the boundary between 
the equal-pairing core and polarized normal domain. It 
should be emphasized that the satellite peak at zero de- 
tuning is indirect evidence for the FFLO pairing. This 
prediction can be experimentally checked by carefully ex- 
amining the RF spectroscopy in the lower temperature 
region.** Also note that the spectrum presented in Fig. 12 
is in good agreement with that obtained from nonlinear 
response theory. 

5. Concluding Remarks 

In this paper, we have theoretically studied the sta- 
ble superfluid state in strongly interacting trapped Fermi 
systems with population imbalance, based on the the 
single-channel Hamiltonian. We have numerically solved 
the BdG equation coupled with the regularized gap equa- 
tion and the number equation in the BCS-BEC crossover 
regime under imbalanced spin densities, where the com- 
putation for the higher-energy contribution was supple- 



mented by the LDA. 

The main results are twofold: (1) First, in § 3, we have 
discussed the ground state in the crossover regime under 
population imbalance and presented the quantum phase 
diagram in the l/kpa-P plane. In the weak-coupling 
regime (l/fcpa < 0), it has been found that the result- 
ing pairing field at T = exhibits the FFLO oscillation 
around the edge of the minority component. In partic- 
ular, the pairing field exhibits the oscillation in the en- 
tire region of the system when P approaches Pc- This 
novel pairing state is refiected in the density and local 
magnetization profiles, such as the bimodal structure for 
the minority species. In contrast, the FFLO oscillation 
disappears for all P in the BEC regime, where the re- 
sulting ground state yields the phase separation between 
the balanced pairing domain and the fully polarized spin 
domain. We have found that the spatial variation of the 
pairing field affects the density. For instance, the presence 
of the FFLO-modulated pairing field leads to a partially 
polarized spin density, while the PS state is refiected in 
the fully polarized spins. It has also been shown that 
the FFLO modulation survives even in the unitary limit 
as the proximity effect. We have confirmed that these 
tendencies of the ground-state structure are unchanged 
in a system with a realistic particle number A^~C'(10^), 
which is comparable with recent experiments.^"^ Our cal- 
culations reproduce a PS-like profile in the local magne- 
tization, while the pairing field yields the FFLO modu- 
lation even in the vicinity of the resonance. The period- 
icity and intensity of the modulation increase as l/kpa 
approaches the weak coupling BCS regime. In particu- 
lar, we have found that the periodicity L of the FFLO 
oscillation is well scaled with the coherence length as 
L ^ 3^0, throughout the extensive range of fcra < and 

P<Pc- 

(2) The second part of the present paper has been 
devoted to another observable quantity, the RF spec- 
troscopy. By numerically solving the tunneling current 
derived from linear response theory, the contributions of 
the pairing field to the spectrum have been discussed. 
The clear difference in the resonance shape between the 
BCS and BEC sides reveals the different superfiuid state, 
which can be checked by further experiments. 

Note that in the phase contrast imaging of the local 
magnetization by the MIT experiment,'^ a PS-like profile 
was observed at the unitary limit l/fci,'a = 0. Our results 
presented in Figs. 6 and 8 are in good agreement with this 
profile, which implies that the pairing field exhibits the 
FFLO oscillation at the edge of the cloud. Also, we have 
demonstrated that as the system approaches the weak 
coupling BCS regime, the FFLO modulation covers the 
entire region of the system, particularly at P ^ Pc- In 
a previous experiment,^ the quantum phase transition 
in the BCS side of the resonance, l/fcpa ~ —0.4, has 
already been observed, which is a favorable condition for 
the detection of the FFLO. Further detailed analysis may 
catch the signature of the FFLO pairing via the density 
and RF spectroscopy results, as we have described in the 
present paper. 

Finally, we comment on the thermodynamic stability 
of the FFLO phase against the increase in T. In Fig. 13, 
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Fig. 13. (Color online) Typical T-P phase diagram in the weak 
coupling BCS side of a resonance (l/fepa = — 0.75). The dashed- 
dotted line denotes Tc for the non-oscillating BCS state. Empty 
circle is the Lifshitz point. 



we present a phase diagram in the T-P plane in the BCS 
side. The phase diagram yields multiple superfluid phases 
composed of the FFLO state and the nonoscillating BCS 
state. All three phase transition lines between the BCS, 
FFLO, and normal states are of the second order, and 
these lines meet at the so-called Lifshitz (L) point. ^"^ The 
BCS-FFLO hne starts from P = r = 0, implying that the 
ground state is always the FFLO state when P ^ and 
T is low. The BCS state only appears at higher values 
of T. This second-order phase transition via the FFLO 
state may be realized in the weak-coupling regime, while 
in the strong-coupling regime the first-order transition 
is predicted by a LDA calculation^^ in which the Lif- 
shitz point is replaced by the tricritical point. We also 
note that this L point exhibits the universal temperature 
Tl/Tco ~ 0.6 in the weak-coupling regime, independent 
of the coupling constant l/fcpa. Hence, all three second- 
order lines are uniquely determined with a fixed L point. 
A similar phase diagram has been proposed even in the 
absence of the trap potential. 

One of the main outcomes in the current work is that 
the FFLO-modulated pairing field survives in the whole 
region in the system approaching the weak-coupling limit 
(l/fcpa — > — cxd). The calculations presented here have 
been performed in a system restricted to cylindrical ge- 
ometry. The FFLO oscillating pattern in a fully three- 
dimensional system without any restriction, such as an 
elongated cigar-shaped or disk-shaped trap, is still open 
to question, and should be further explored in future. 
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Appendix A: Regularized BdG Equation 

We start with the original Hamiltonian in eq. (2). Here, 
it is convenient to introduce a spinor in the Nambu space. 



*(r) = [V;T(r),^ (r)] 



(A-1) 



Applying the standard mean-field approximation to the 
interaction part of the above Hamiltonian, the effective 
Hamiltonian can be derived as 



Hmy — f + 
where 



dr J dr'*t(r)/C(r,r')*(r'), 



/CT(r,r') 
A*(r,r') 



A(r,r') 
-q(r,r') 



(A-2) 
(A-3a) 

/C,(r,r') = Hi°^S{r - r') + >V_,(r, r'). (A-Sb) 

The mean-field quantities, the pairing potential A(r,r') 
and the Hartree potential Wo-lr, r'), are defined by 

A(r,r') = C/(f)(V'x(r')V'T(r)), (A.4a) 

Wa{r,r') = U{f)p„{T) = C/(f)(^Ai(r)^,(r)). (A-4b) 

is the c-number including the condensation and 
Hartree energies and f = |f| ~ |r — r'| is the relative 
coordinate. 

The Bogoliubov transformation of the spinor ^ (r) into 
the quasi-particle basis r]^ = [T]„,/i ,f]l J"^ is defined by 



*(r) 



E 



v^{r) <(r) 



Here, the creation and annihilation operators of the 
quasi-particles, rjl ^ and rj^^, obey the fermionic com- 
mutation relations. The quasi-particle wave function in 
the matrix form u^{r) satisfies the orthonormal condi- 
tion, 

J drul{r)u^'{r) S^.y. (A-6) 

Also, we write the completeness in a matrix form, 

J2Mr)ulir')^S{r-r'). (A-?) 

Then, we assume that the mean- field Hamiltonian can 
be transformed into the diagonalized form, TYmf = + 
YlivYlia^^"^ '^'l,a"']v,a- To this end, we obtain the BdG 
equation 



dr'^(r,r')u^(r') = u^(r) 



'■■V 









(A-8) 



The above BdG matrix /C yields double eigenstates for 
hyperfine spins. To see this, we set the eigenfunction with 
the up spin as ^p'p = \u,y^Vi^ ^ with an eigenvalue of 
e^v^: t^p]}^ ^ei)^^p]}^ . It is found that the BdG equation 
(A-8) simultaneously has eigenstates for down spins of 



ip];' = [—?;*, u*] with the eigenvalue —e 
one can solve the BdG equation 



(i) 



. Therefore, 



(A-9) 



using the eigenfunction (p^, = [m^, d^]^ and the eigenstates 
corresponding to those for up and down spins. How- 
ever, we should emphasize that e^^ ^ si^^ because the 
finite mismatch of the Fermi surface causes the breaking 
of the time-reversal symmetry, /C(r, r') ^ — f2/C*(r, r')f2, 
where Tj is the jth Pauli matrix. 
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Using the eigenstates in cq. (A-9), the pairing field 
defined in eq. (A-4a) can be explicitly expressed as 



A(r,r') = U{f)Y, [uArK{r')fP - <(r)n.(r')/(^) 



C/(f)^w.(r)<(r')/., 



(A 



^^^) = {Vu',ivli)Suy, and 
The summation in the gap 



for the Fermi distribution functions /i^' = /(ei^^) = 

ivl ^Vi^^)Suy, f^^^ = /(-£; 
/.^ /(£;.) = 1/(6^^/^ + 1 
equation (A- 10) is performed for all the eigenstates with 
both positive and negative eigenenergics. Then, the par- 
ticle density in each spin state is given from Eqs. (A-4b) 
and (A-5), and is described in cq. (9). The BdG equation 
(A-9) is now self-consistcntly solved under the mean- field 
conditions, given by Eqs. (A- 10) and (9), for the fixed to- 
tal particle number: N = J2a- — I 

Now, let us derive the explicit expression for an equal 
mixture of two-component fermions, i.e., /C| under 
(5/i = 0. Then, the BdG matrix (A-3a) yields the time- 
reversal symmetry, /C(r,r') = — f2/C*(r, r')f2. It hence 
follows that a positive eigenvalue i?^ having the eigen- 
function ip^ = [uj,, Vi,]'^ in eq. (A-9) is always accompa- 
nied by the negative eigenvalue —E^, with the eigenfunc- 
tion —iT2<pl = [— u*,zt*]"^. Then, the mean-field quanti- 
ties shown in Eqs. (A- 10) and (9) can be reduced to the 
standard form, 

A(r, r') = U{f) ^ «.(r)<(r')[2/. - 1], (A-lla) 

E^>0 

[K(r)|V. + k(r)|'(l -/.)]. (A-llb) 



where JFrog is the regular part of the anomalous average 
^(R,f) = (V'i(r')V'TW) (R=(r + r')/2 is the center-of- 
mass coordinate). One way to remove the divergent term 
is to replace the original contact interaction g with the 
pseudopotential''^^^^ and then, the formal expression for 
j^Q^he pairing field is given by 

A(r) = lim ^ [f (V'^(r')VT(r))] - (A-14) 

To give a straightforward expression for the regular- 
ization operator lim^^o c^f ['Ji we introduce the single- 
particle Green's function, G^(r, r') = {r\H^^\r') with 
6)1 — 0. This function yields the same nature of divergence 
as A in the limit of f ^ 0: G^(R, r) = ^ + G]^s(R) + 
0{'r). Here, the divergent contribution of the anomalous 
average in the right-hand side of eq. (A-14) is canceled 
out by the irregular part of the single-particle Green's 
function, which allows one to introduce an explicit en- 
ergy cutoff. By employing the LDA, we finally obtain 
the regularized gap equation, 



E 



Hereafter, we use the definition in the case of a system 
without the time-reversal symmetry. 

At low temperatures, the collisions between atoms can 
be described in terms of s-wave scattering, U{r)= g6(r). 
Here, g = 4'!T'^a/M is the coupling constant, and in gen- 
eral the interaction can be expressed using the dimen- 
sionless parameter k-pa, where fcp = ^/2ME-p is the Fermi 
wave vector of a noninteracting Fermi gas with estimated 
Fermi energy Ey- Then, the BdG equation (A-9) is re- 
duced to the following local form: 



/C^r) 
A*(r) 



A(r) 

-/C|(r) 





M,y(r) 


= E, 











(A-12) 



where A(r, r') ==(5(f)A(r) and /C<,(r, r') = (5(f )/C^(r). 

Here, there are two singular contributions to the BdG 
equation (A-12) and gap equation (A- 10). First, the 
Hartree potential gp^ diverges at the unitary limit, which 
is neglected throughout this paper (sec the text in § 2.1, 
for dcrtails) . The second singular behavior arises from the 
fact that the contact interaction leads to the UV diver- 
gence of the pairing field defined in eq. (A- 10), where 
the leading term of A(r, r') behaves as — MA(r)/47rf at 
f ^ 0: 



A(r,r') 



47rr ^ ' 



-(?J-,eg(R,f)+0(f), (A- 13) 



A(r)=5(r)^«,(r)<(r)/„ 



(A- 15) 



where the renormalized coupling constant ^(r) is given 

by 



1 _ 1 ^ Mfce(r) 
fl(r) g 27r2 



1-^1-^4^^ (A-16) 
2A;c(r) fcc(r) - fcpCr) J ^ ' 

Here, fcF(r) and kc{r) are the local wave vectors defined 
by the local Fermi and cutoff energies, respectively (see 
eq. (7)). 

Appendix B: Matsubara Green's functions 

The Matsubara Green's function in the Nambu space 
is defined using the imaginary time r as 

^(1,2) = -(t, [*(l)*t(2)]y (B-1) 

Here, we introduce the notation 1 = (ri, ri). The Green's 
function in a 2 X 2 matrix form obeys the Gor'kov equa- 
tion 



^F(r) 



88 



/ 



d2 



d_ 



fo<5(l,2)-K;(l,2) 



g(2,l')=fo(5(l,l'),(B-2) 



for a 2 X 2 unit matrix fo- Also, we set ^(1, 2) = S{ti — 
T2)Ai(ri, r2). From the BdG equation (A-8), ^(ri,r2) 
may be expanded as 



^(ri,r2) = ^u^(ri) 



-(T) 








uUv2). (B-3) 



The Green's function becomes diagonal in the repre- 
sentation that /C is diagonal. Hence, ^(1,2) may be ex- 
panded as 

1 --"(-i-=)^ti.(ri)4(ia;„)ut(r2). 



5(1,2) = - > e 



(B-4) 



Its Fourier component G^^^\e) is obtained easily as 
[iun + ei^)) 



(B-5) 
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It is convenient to introduce the Fourier 
transform with respect to r, Q{1,2) = 
/3^^ e~"^"('^i~'^^)^(rir2, iwn), whose coefficient 
is obtained from cq. (B-4). From the eigenfunctions 
and energy of the BdG equation, {Ui,,Vi,) and E^, 
^(rir2,iw„) is given by 



yil{rir2, ILOn) = } —. 



M t/(ri)M*(r2) 
iu)n - El, 



r I • ^ ^^i'(riK(r2) 



(B-6a) 
(B-6b) 
(B-6c) 

(B-6d) 



where Qij denotes the element of the 2x2 matrix 
Q. 
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